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Abstract 


Particle dispersion and the influence particle momentum exchange on the properties of a turbulent 
carrier flow in micro-gravity environments challenges understanding and predictive schemes. The 
objective of this effort has been to develop and assess high-fidelity simulation tools for predictmg 
particle transport within micro-gravity environments suspended in turbulent flows. The compu- 
tational technique is based on Direct Numerical Simulation (DNS) of the incompressible Navier- 
Stokes equations. The particular focus of the present work is on the class of dilute flows in which 
particle volume fractions and inter-particle collisions are negligible. Particle motion is assumed 
to be governed by drag with particle relaxation times ranging from the Kolmogorov scale to the 
Eulerian timescale of the turbulence and particle mass loadings up to one. The velocity field was 
made statistically stationary by forcing the low wavenumbers of the flow. The calculations were 
performed using 96 3 collocation points and the Taylor-scale Reynolds number for the stationary 
flow was 62. The effect of particles on the turbulence was included in the Navier-Stokes equations 
using the point-force approximation in which 96 3 particles were used in the calculations. DNS 
results show that particles increasingly dissipate fluid kinetic energy with increased loading, with 
the reduction in kinetic energy being relatively independent of the particle relaxation time. Viscous 
dissipation in the fluid decreases with increased loading and is larger for particles with smaller re- 
laxation times. Fluid energy spectra show that there is a non-uniform distortion of the turbulence 
with a relative increase in small-scale energy. The non-uniform distortion significantly affects the 
transport of the dissipation rate, with the production and destruction of dissipation exhibiting com- 
pletely different behaviours. The spectrum of the fluid-particle energy exchange rate shows that 
the fluid drags particles at low wavenumbers while the converse is true at high wavenumbers for 
small particles. A spectral analysis shows that the increase of the high wavenumber portion of the 
fluid energy spectrum can be attributed to transfer of the fluid-particle covariance by the fluid tur- 
bulence. This in turn explains the relative increase of small-scale energy caused by small particles 
observed in the present simulations as well as those of Squires & Eaton (1990) and Elghobashi & 
Truesdell (1993). 
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1 Introduction 


The interaction of solid particles or liquid droplets with gas-phase turbulent flows controls the 
performance of many engineering devices and is important in natural processes as well. Exam- 
ples include the combustion of pulverised coal or liquid sprays, transport of particulate solids, 
gas-phase reactions controlled by particulate catalysts, dust storms, and atmospheric dispersal of 
pollutants. Especially relevant to the present effort are particle-fluid interactions in micro-gravity 
environments. In each of these areas an increased understanding of the fundamental phenomena 
that drive the complex interactions between the particle cloud and turbulent earner flow is needed 
to ultimately improve the design of engineering devices in which these flows occur. 

Within the vast array of applications encompassed by two-phase flows, the particular interest 
of the present work is on the interaction with a turbulent carrier flow of a dilute dispersed phase 
of particles with densities substantially larger than the carrier flow, i.e., P 2 I Pi ~ C^IO 3 ) where 
P 2 and p\ denote particle and fluid density, respectively. For dilute flows the volume fraction of 
the particles, a 2 > is small. However, the particle mass loading, <f> = 0 . 2 P 2 .I P\ , can be large enough 
such that momentum exchange between particles and fluid results in a significant modulation of 
the turbulence, typically referred to as two-way coupling (e.g., see Crowe et al. 1996). 

For this class of two-phase flows the modulation of turbulence by particles is complex and still 
not well understood. For example, experimental measurements in shear flows, e.g., particle-laden 
jets and boundary layers, have shown that turbulence velocity fluctuations may be either increased 
or decreased due to the modulation of the flow by heavy particles (e.g., see Tsuji et al. 1984, 
Modarress et al. 1984, Shuen et al. 1985, Fleckhaus et al. 1987, Hardalupas et al. 1989, Gore & 
Crowe 1989, Rogers & Eaton 1991, Kulick et al. 1994). In turbulent shear flows it is often difficult 
to separate the direct modulation of the turbulence due to momentum exchange with particles 
from the indirect changes occurring through modification of turbulence production mechanisms 
via interactions with mean gradients. Furthermore, it is usually difficult in experiments to isolate 
the effects of different parameters on measurements. 

Numerical simulation offers another approach for examining the interactions between particles 
and turbulence and the modulation of turbulence by particles. For particle-laden flows traditional 
approaches relying on solution of the Reynolds-averaged Navier-Stokes equations require empir- 
ical input, principally the prescription of turbulence properties along particle trajectories. For 
applications of one-way coupling, i.e., no modulation of the flow, various modeling approaches 
have been developed which adequately describe dispersion in simple flows, though the central 
problem of prescribing Lagrangian quantities along particle trajectories remains an open ques- 
tion (e.g., see Simonin et al. 1993, 1995). For applications of turbulence modulation by particles 
turbulence quantities such as the kinetic energy and dissipation rate are modified directly by the 
particles (e.g., through the appearance of source terms in the transport equations) as well as indi- 
rectly through changes which particles causes in turbulence dynamics (e.g., see Squires & Eaton 
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1994). Thus, the empirical input required for Reynolds-averaged approaches at the present time 
makes it difficult to obtain a fundamental understanding of turbulence modulation. 

The most sophisticated numerical approach for examining particle-turbulence interactions is 
direct numerical simulation (DNS). In DNS the Navier-Stokes equations are solved without re- 
sorting to ad hoc modeling at any scale of motion. The primary advantage for calculation of 
particle-laden flows is that turbulence properties along particle trajectories are directly available. 
For applications of one-way coupling several studies have examined heavy particle transport in 
isotropic turbulence (e.g., see Deutsch & Simonin 1991, Squires & Eaton 1991a, Elghobashi & 
Truesdell 1993, Wang & Maxey 1993). 

For two-way coupling DNS has been applied to particle-laden isotropic turbulence and has 
demonstrated that the distortion of the turbulence is not uniform and is dependent upon the par- 
ticle relaxation time (e.g., see Squires & Eaton 1990, Elghobashi & Truesdell 1993). Squires & 
Eaton (1990) considered particle motion in the Stokes regime in which gravitational settling was 
neglected. Computations were performed using both 32 3 and 64 3 grids at Taylor-scale Reynolds 
numbers of 35 in which a steady, non-uniform body force was added to the governing equations 
in order to achieve a statistically stationary flow. Particle sample sizes up to 10 6 were used m the 
simulations. Mass loadings from zero (one-way coupling) to unity were considered for a series 
of particle relaxation times varying from 0.3r fc to llr fc where r fc is the Kolmogorov time scale. 
For a Stokes drag law without gravitational settling it is straightforward to show that particles will 
globally dissipate turbulence energy. Squires & Eaton (1990) found that the overall reduction in 
turbulence kinetic energy for increasing mass loading was insensitive to the particle relaxation 
time. They also showed a strong preferential concentration of particles into regions of low vortic- 
ity and/or high strain rate (see also Wang & Maxey 1993). For cases of turbulence modulation. 
Squires & Eaton (1994) attributed the non-uniform distortion of the turbulence energy spectrum 
by particles to preferential concentration. 

Elghobashi & Truesdell (1993) examined turbulence modulation by particles in decaying isotropic 
turbulence using resolutions of 96 3 for the Navier-Stokes equations and 34 3 particles. Particle mo- 
tion in Elghobashi & Truesdell (1993) was governed by the equation derived by Maxey & Riley 
(1983). They found that for the large density ratios considered in their simulations particle motion 
was influenced mostly by drag and gravity. Elghobashi & Truesdell (1993) found that the coupling 
between particles and fluid resulted in an increase in small-scale energy. The relative increase in 
the energy of the high wavenumber components of the velocity field resulted in a larger turbulence 
dissipation rate. They also found that the effect of gravity resulted in an anisotropic modulation 
of the turbulence and an enhancement of turbulence energy levels in the direction aligned with 
gravity. Furthermore, in the directions orthogonal to the gravity vector a reverse cascade of energy 
from small to large scales was observed. 

While the work of Squires & Eaton (1990) and Elghobashi & Truesdell (1993) have advanced 
our understanding, the effect of turbulence modulation by particles is not frilly resolved. For ex- 


2 


ample, in the transport equation for turbulence kinetic energy the coupling between particles and 
turbulence yields an additional term which accounts for the energy transfer imparted from the parti- 
cles to the fluid. Both Squires & Eaton (1990) and Elghobashi & Truesdell (1993) have shown that 
the distortion of the turbulence energy spectrum is sensitive to quantities such as the particle relax- 
ation time. This implies that the energy transfer from particles to turbulence acts non-uniformly 
across the spectrum. However, its overall behaviour, nor its spectral distribution, is available from 
previous investigations. Both the global value as well as spectral distribution are important not only 
for increasing fundamental understanding but also development of engineering models. Relevant 
in this regard is the work of Baw & Peskin (1971) who have previously considered the spectral 
modulation of turbulence by particles. They showed that the effect of particles is to decrease the 
energy at high wavenumbers more than that at low wavenumbers. Their analysis, however, contra- 
dicts the results of both Squires & Eaton (1990) and Elghobashi & Truesdell (1993). 

The objectives of the present work are to investigate turbulence modulation by particles within 
micro-gravity environments. The particular flow field of interest being isotropic turbulence. In 
isotropic turbulence there is no production and therefore from a given initial condition, the flow 
decays over time, i.e., turbulence time and length scales increase. In decaying turbulence the evo- 
lution of quantities such as the kinetic energy and dissipation rate exhibit a dependence on initial 
conditions. The decay of the flow and dependence on initial conditions complicates interpretation 
and analysis of both particle motion and turbulence modulation. An alternative to simulation of 
decaying turbulence is calculation of flows made statistically stationary through an addition of a 
body force to the Navier-Stokes equations in which the force is added to the low wavenumber com- 
ponents of the velocity field. Statistically stationary flows can be advanced to an equilibrium in 
which particle motion, and the effect of particles on the flow, are independent of initial conditions. 
Time and length scale ratios of the turbulence relative to the particles are also stationary. The forc- 
ing scheme used in this work is that developed by Eswaran & Pope (1988), who have shown that 
the small scales of the velocity field are insensitive to the energy input from the forcing. Discussion 
of computation of two-way coupling in DNS is presented in §2. The point-force approximation is 
used in this work to account for momentum transfer between particles and turbulence and impor- 
tant issues relevant to this approach are discussed. An overview of the simulations is also presented 
in §2 with evolution of statistical quantities and the spectral analysis in §3. A summary of the work 
may be found in §4. 
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2 Simulation Overview 

2.1 True direct numerical simulation of two-phase flows 

True direct numerical simulations of fluid flows loaded with heavy particles require that one re- 
solves the standard Navier-Stokes equations for the fluid: 

^1 = 0 (0 

dxi 


du u dui t i _ 1 dpi d 2 «i,j 

~dt + Ul ' k dx k ~ Pi dxi 1 dx k dx k ' 


( 2 ) 


The effect of the particles on the fluid is formally taken into account through the boundary condi- 
tions on the surface of each particle, 


Ul i ( X) t) = w?(x, t ) for all points on the particle surface, i.e., (x, t) € , (3) 

where < is the instantaneous velocity of the surface of particle n (the n superscnpt is used 
throughout this work to denote properties of a single particle). For rigid particles m translation, 
w n is the same everywhere on O n and is equal to the velocity at the centre of the particle. The 
subscripts 1 and 2 denote the fluid and particle phases, respectively. Thus, in (1) and (2), u ui is 
the i th component of the fluid velocity, p l the fluid pressure, and p x and v x the fluid density and 

kinematic viscosity, respectively. 

Simultaneous to the solution of (l)-(3), particle trajectories are computed using Lagrangian 
tracking, the force acting on the particle being computed by direct integration of the simulated 
fluid stress on the particle surface: 


dvSi f 

P2—tt = mi + 
at J n n 


—p\5ij + v\ 


du i,i 
dx j 


rij dto . 


(5) 


The displacement of particle n is x^, the particle density is p 2 . The outward pointing normal 
to the surface fi" is n jt and g t is the acceleration of gravity. The approach outlined in (l)-(5) 
requires that the velocity field around each particle be accurately resolved and is viable only for 
calculations with 0(10) particles (e.g., see Unverdi & Tryggvason 1992). In a turbulent flow with 
a large ensemble of particles having diameters on the order of the Kolmogorov length scale this 
approach is not feasible owing to its enormous computational cost. Thus, approximations are 
required in order to account for the effect of particle momentum exchange on the flow. 
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2.2 Point-force approximation 

Saffinan (1973) showed that the perturbation in the fluid due to the presence of a particle decays 
as the sum of two contributions, one as 1/r (long-range) and the other as 1/r 3 (short-range). For 
particles small relative to the smallest length scales of the flow, and for particles separated by 
a distance L large compared to their diameter d, the most important interactions are long-range 
(e.g., see Koch 1990). Neglecting short-range interactions, e.g., particle wakes, is justifiable for 
particles with diameters smaller than the Kolmogorov length scale of the flow field undisturbed by 
the presence of the particle since in that case short-range perturbations are dissipated by viscosity. 
The focus of this work will be on larger scales in which long-range interactions are dominant. 
The Navier-Stokes equations can then be written for the fluid everywhere in the domain with the 
influence of particles taken into account by Dirac distributions of the force, /£, applied to the fluid 

by each particle, 


where 




05 

H ^ 

II 

O 

(6) 

dttl.i 

dt 

du i,, 

' + U\ k a 

ox k 

1 dpi , d 2 u ui t 1 f 

■ = — i^l a o + „ Jc,z i 

p x dxi ox k ox k Pi 

(7) 


- fa 

pi 

= /d,t( a '2,i)^( X * — X 2,i) • 

Pi 

(8) 


Hie IOrce IS opposite oi uiai appuvu ivj 

and Maxey & Riley (1983), the forces acting on a particle can be considered to arise from three 
contributions. The first contribution, f 0 , represents the virtual force that would apply on a fluid ele- 
ment that coincides with the particle position, i.e, pressure forces and viscous stresses. The second 
contribution, f 6 , arises from the perturbation of the fluid flow due to the presence of the particle. 
For a rigid sphere of diameter d in translation, this perturbation of the surrounding unsteady non- 
uniform flow results in the drag, added mass, and Basset history forces. The third contribution is 
gravitational settling. For spheres with density p 2 large compared to the fluid density p u f b > f a 
and reduces to the drag force. The particle equation of motion can then be written for a single 

particle n as 

i n 

(9) 


P 2 


du^ ; 3 

’ = ~fd,i + Pi9i — _ 7Pi 


dxli_ n 

dt “ 2 ’‘ ’ 


v" I + P29z = -P2 


T P 


+ p29i ■ (1®) 


dt Ja ’ 1 ' r ‘" y ' 4 ri d 
The local drag coeflicient in (10) is C£ and may be expressed in terms of the particle Reynolds 
number Re™ as (Clift et al. 1978) 


Ck 


24 

Re™ 


1 + 0.15i?e” 


0.687 


Re n P = 




< 800 , 


( 11 ) 


5 



( 12 ) 


where r p is the particle relaxation time, or time constant, defined as 


T„ = - 


4 P2Cl.n1 

3 pi d 1 rl 


d 2 P 2 _2 

I 81 / 1 P 1 1 + 0.15Ref 687 


The local instantaneous relative velocity between particle n and the surrounding fluid is = 
u^i - v-ii, where is the fluid velocity at the position of particle n of the flow field locally 
undisturbed by the presence of the particle (Gatignol 1983, Maxey & Riley 1983). The expression 
for the drag force as written above is applicable to particles having diameters smaller than the 
Kolmogorov length scale, i.e., d < rj. If d is comparable to 77 , Faxen terms should be taken 
into account (Gatignol 1983, Maxey & Riley 1983). Gravitational settling results in the crossing 
trajectories effect which strongly influences particle dispersion (Csanady 1963, Squires & Eaton 
1991b). Elghobashi & Truesdell (1993) have shown that in two-way coupling gravity complicates 
turbulence modulation. In the current study the parameters varied are particle relaxation time and 

mass loading, the effect of gravity is not examined. 

For a flow containing N p particles, the fluid velocity required in (10) is that locally undis- 
turbed by the presence of particle n, but taking into account the disturbances created by all other 

_ i) particles in the flow. This in turn requires that to determine the motion of each particle, a 
total of N p flow fields are required. Only in the limit of one-way coupling is u”, identical to that 
in a single-phase turbulent flow. 

Durlofsky et al. (1987) showed that the exact resolution of two-phase flows two-way coupled 
with N p particles under the point-force approximation requires inversion of a system of N p non- 
linear equations, which becomes impossible for large N P . The difficulty to simulate, under the 
point-force approximation, turbulence modulation by a large number of particles can be illustrated 
by first considering a flow initially at rest with a single particle that influences the fluid. The 
locally undisturbed fluid velocity is that without the particle. The drag force is then —u 2 t i/r p while 
the fluid velocity resulting from the disturbance created by the particle is the superposition of the 
background flow and the perturbation due to the particle that is, for a fluid at rest containing a 
single particle, a Stokeslet (e.g., see Gatignol 1983). For a flow modified by a large number of 
particles the motion of particle n requires knowledge of the velocity that is locally undisturbed by 
its presence, but takes into account the disturbance created by all other particles in the flow. In 
this case the locally undisturbed fluid velocity for particle n is again that in which the particle is 
not present in the flow but its influence must be maintained in determining the motion of all other 

particles. 

Thus, numerical simulation of two-way coupling requires that there must be as many undis- 
turbed fluid fields as particles. Durlofsky et al. (1987) showed that this requires inversion of a 
system of N p non-linear equations, which becomes impossible for large N p . For practical calcula- 
tions a unique fluid velocity field is desired which is influenced by the entire ensemble of particles 
but may also be used to obtain the velocity and position of each particle, i.e., a velocity field that 
can be considered as the locally undisturbed fluid velocity for obtaining the motion of each particle. 
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2.3 Effective direct numerical simulation of two-phase flows 

One approach to obtaining the locally undisturbed fluid velocities would be to resolve the 
flow field U\ i influenced by the entire ensemble of particles and then subtract the local perturba- 
tion induced by the particle presence. Considering Stokes flows, i.e., particle Reynolds number 
Re v < 1 , in the dilute regime, i.e., inter-particle separations large with respect to the particle di- 
ameter, Saffman (1973) showed that u M is the sum of and the local perturbation, the so-called 
Stokeslet, 


u l,i = ^l.t + 


3 d ' 

4 r . 


v*, + v 


r,J r 2 


with r = | x — X 2 


r = x — x 9 


(13) 


To obtain the locally undisturbed velocity field for each particle, u ” an iteration procedure could 
be developed using (13). For a flow with a large number of particles, however, the computational 
cost becomes prohibitively large. 

For practical purposes with a large sample of particles N p , and for intermediate particle Reynolds 
numbers, it is necessary to assume that for each particle, the locally undisturbed fluid velocity field 
Uii can be approximated by Therefore, the coupling force may be expressed as, 


fn 

I d,i 


P 2 


x 2,i 


— U 


l,i 


U 


T 

P 


P 2 - 


2,i - «l.t 


(14) 


An estimate of the error made in simulations in which the approximation (14) is employed can be 
obtained using (13). Defining A u^ t as the error, (13) shows that 




Ui,i = 


3d 

4r 


v ■ + v • — - - 

u r, 1 ' u r,J r 2 


(15) 


In actual computations, the distance r between particle n and the grid nodes (where the locally 
undisturbed fluid velocity, is approximated by that from the DNS, u hi ) is of the order of 
the mesh size. In simulations, the distance r between particle n and the grid nodes is of the 
order of the mesh size Ax. Thus, on average, the relative error resulting from the use of (14) is 
0(d/Ax) and, in addition to the restriction d-Cr/ imposed by the point-force approximation and 
neglect of the Faxen contributions in the particle equation of motion, the condition imposed by 
the approximation of the locally undisturbed velocity u” j by U\a requires that d Ax. These 
constraints are compatible since in a DNS calculation, 77 is of the order of Ax. This constraint is 
consistent with the assumption of the point- force approximation which retains only the long-range 
fluid-dynamic interaction. Since in Direct Numerical Simulation (DNS) 77 is of the order of Ax, 
the above constraint is consistent with the assumption of the point-force approximation to retain 
only the long-range fluid-dynamic interaction. 

Furthermore, the identification between and can be understood as the negligence of the 
local perturbation of particle n on itself with respect to the influence of all the other particles on 
particle n. If the N p particles are uniformly distributed, separated, on average, by a distance L, 
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they form a regular “particle-box”. If one considers one particle n and looks at the perturbations 
created by its nearest particles (6 at a distance L, 12 at a distance v/2T and 8 at a distance \/3 L), the 
resulting perturbation dues to these particles on particle n writes as a sum of Stokeslets, assuming 
that these Stokeslets can be based on the same u\ t i, u”, and u^- The factor is. 


3d 

4 L 


6 + 


12 _ 8 _ 
V2 + y/Z. 


(16) 


Therefore, the Stokeslet due to particle n is negligible compared to those created by the nearest 


particles when: 


d 3d 

-r- « 77 

Ax 4 L 




15 d 

~T 


(17) 


Thus including or not the effect of the Stokeslet due to particle n on itself is negligible when 


L <C 15Ax 


(18) 


2.4 Direct numerical simulation of isotropic turbulence 

The method used to obtain the fluid velocity is based on the direct numerical simulation technique 
developed by Rogallo (1981) in which dependent variables are expanded in Fourier series and the 
flow is represented in a cubic domain of volume Lf )0X = (27t) 3 with periodic boundary conditions. 
Exact integration of the viscous terms is performed using an integrating factor and the non-linear 
terms are calculated in physical space. The discretised equations are time advanced using a second- 
order Runge-Kutta scheme. The reader is referred to Rogallo (1981) for further details on the 
method. 

The method is applied to computation of homogeneous isotropic turbulence. Isotropic turbu- 
lence is non-stationary since in the absence of a production mechanism turbulence decays. As 
discussed in § 1, the continual evolution of turbulence quantities complicates analysis and interpre- 
tation of decaying turbulence. Lack of a statistically stationary flow in particle-laden turbulence is 
even more complex since the ratio of the particle relaxation time to fluid time scales changes as the 
flow evolves. To alleviate these complications, a spatially non-uniform, time-dependent body force 
(or acceleration) was added to the low wavenumber components of the velocity field to maintain a 
statistically stationary flow. 

The large scales were forced using the scheme developed by Eswaran & Pope (1988) and 
Yeung & Pope (1989) and is based on an Uhlenbeck-Omstein stochastic process that determines 
an acceleration for each of the three velocity components and for each non nul wavenumber mode 
within a shell in spectral space of radius Kp. The complex-valued body force is added to the 
momentum equations at each time step. The evolution of one component of the body force / at 
time level n + 1 is obtained via 

fn+i = /n(l - dt/T F ) + e^/2a 2 dt/T F , (19) 
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where e is a random number taken from a Gaussian distribution of zero mean and unit variance. 
The characteristic time of the process, Tp, as well as the forcing variance, c, were equal to 1.4 
and 0.033, respectively. A simple modification of the forcing was adopted in this study since 
the Uhlenbeck-Omstein process assumes that the time step dt remains small with respect to the 
characteristic time of the energetic turbulent structures, i.e., with respect to the Eulerian time scale 
T e (defined using the Eulerian fluid velocity and longitudinal integral length scale, L e ). It is also 
assumed that the time step dt remains large with respect to the small scale motion, characterised 
by the Kolmogorov time scale r^, which ensures that the forcing is independent of small scale 
motions (Eswaran & Pope 1988). Therefore, the forcing acceleration was maintained constant 
over 2 t£ with a ratio T e /T^ ~ 10. This small modification of the forcing scheme requires longer 
averaging periods in order to ensure isotropy of the turbulence. 

For the simulations reported in this paper the turbulence was resolved using 96 3 collocation 
points with a forcing radius Kp equal to eight. This resolution provides an adequate separation 
between the forced modes and small scales and a large enough Reynolds number to obtain a sepa- 
ration between the peak of the energy and dissipation spectra (Yeung & Pope 1989, Boivin 1996). 
The maximum wavenumber of the simulation k max is \/2/3N with N the number of grid points 
in each direction (Rogallo 1981). As shown by several investigators values of k max T] greater than 
one ensure accurate resolution of small-statistics as well as accurate interpolation of fluid veloc- 
ities (e.g., see Eswaran & Pope 1988, Balachandar & Maxey 1989). The Courant number of the 
computations was 0.5 in order to minimise time stepping errors and ensure accurate resolution of 
the small scales (see Eswaran & Pope 1988 and Yeung & Pope 1988 for further discussion). 

Table 1 summarises the main characteristics of the reference fluid flow, i.e., without influence 
of the particles on the turbulence, that was obtained following a time development required for the 
flow to become independent of its initial conditions. The statistics were accumulated over a period 
of roughly 7r e . In Table 1, q\ and e x are the fluid kinetic energy and dissipation rate, respectively, 
and are used to form the Eulerian time macro-scale, t £ = q\/e i- The Reynolds number, Re\, is 
based on the Taylor micro-scale A, defined as 


9? 


= 2 < 


'2 


>1 n ^ ^1 5 


A = 


‘ 15^i < ti] 2 >i 
Sl 


Re \ = A 


\/< U? >1 


Vi 


( 20 ) 


(the 1 superscript denotes a fluctuating quantity obtained by subtraction from the mean). The kinetic 
energy and dissipation rate are related to the energy spectrum as 


/»oo r°o rcG 

q?= E{k)dk e x = / D{k)dk = / k 2 E(k)dk. (21) 

Jo Jo Jo 

The Lagrangian integral time scale shown in Table 1, t[, is obtained from integration of the La- 
grangian autocorrelation. The Kolmogorov time scale is denoted r* = (^i/ci) 1 '' 2 - Averages of 
turbulence quantities obtained over the computational volume are denoted < • >i. With a ratio 
L) l L bo x around 0. 1 5 , the computational domain contains an adequate sample of energy-containing 
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Table 1 : Flow field parameters 


— 3 — 

i/i q{ 

E\ L C fjL e Q L e f/Lb 0X 

Re x 

kmaxV 

T e 

T £ T\ T X T k 

0.015 7.00 

5.70 1.98 0.148 

62 

1.26 

0.43 

1.23 0.35 0.32 0.051 


eddies to avoid problems due to imposition of periodic boundary conditions (Yeung & Pope 1989). 
Moreover, the ratio of the longitudinal length scale, L), to the transverse length scale, L* g , is close 
to two, in agreement with isotropic relations. Finally, also shown in the table is r A , the time scale 
representative of the Taylor length scale A. It is formally defined with the time scale relation at 
wavenumber k = 1/A valid in the inertial subrange, r A = (A 2 /ei) 1/3 (Hrnze 1975). 

The energy spectrum for the reference flow is plotted in Figure 1 . The energy at low wavenum- 
bers differs from that measured in the grid-turbulence experiment of Comte-Bellot & Corrsin 
(1971) The experimental values are the measurements using a two inch grid at a location where 
Re x = 65, close to Re x = 62 of the reference flow in the DNS. The energy at high wavenum- 
bers follow remarkably well the experimental data, which, considering the comprehensive study 
of Eswaran & Pope (1988), illustrates that the forcing does not adversely affect the small scale 

motion (see Boivin 1996 for further discussion). 

2.5 Numerical implementation of the coupling force 

Properties of the particle cloud were obtained by solving (9) and (10) for a large ensemble of 
particles. A second order Runge Kutta scheme was used for advancement of the particle velocity 
and displacement. Interpolation of fluid velocities to particle positions was performed using third- 
order Lagrange polynomials. Numerical experiments have shown that the scheme is accurate for 
interpolation of quantities such as the fluid velocity (e.g„ see Yeung & Pope 1988, Balachandar & 
Maxey 1989, Boivin 1996). 

Two schemes to incorporate the coupling force in the fluid momentum equations were consid- 
ered. This process is a projection of the coupling force, defined at the particle position, onto the 
grid. The first scheme, the so-called Particle-in-Cell (PIC) method, represents the coupling force 
f c i as proportional to the accumulation offerees /£ (Eqn. 14) induced by each particle n surround- 
ing a node P (with volume V = L box /N) on which the fluid velocity is calculated (e.g., see Crowe 

1982): 

V jP(x,y,z) = Hi , {21) 

n in V 

where a is a constant of proportionality to be defined. In the second scheme, rather than a sum- 
mation of /jv around a node P, the force exerted by each particle on the fluid is projected onto the 
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grid 


Vf?} x ' y ' z) = a ^ P ( x ’ y ' z V ■ (23) 

n in V 

The weights in the projection operation in (23) can be based on the cell volumes as in Squires & 
Eaton (1990) or on the distances between particle n and the eight nearest grid nodes as employed 
by Elghobashi & Truesdell (1993). Both approaches yield similar results (Boivin 1996); therefore, 
only the results in which the weights are based on cell volumes are presented below. 

The two schemes (22) and (23) were tested using an instantaneous fluid velocity field which 
was interpolated to 96 3 random particle locations within the computational box using third order 
Lagrange polynomials. Fluid velocities obtained from the interpolation step were then projected 
back onto the grid using (22) and (23). Figure 2 shows the energy spectra of the fluid velocities 
resulting from this procedure. The second method (23) recovers much more of the kinetic energy 
following the interpolation and projection steps as compared to the PIC scheme, e.g., a decrease in 
the initial kinetic energy of only 2% using (23) compared to a reduction of 38% when using (22). 
The figure also shows that the high wavenumber end of the energy spectrum is more accurately 
recovered using (23) compared to the PIC scheme. Also shown in Figure 2 is the energy spec- 
trum of the fluid velocities following interpolation of the initial field onto a mesh shifted by half 
a grid cell (dotted curve in the figure). This curve shows that there is a relatively small filtering 
of high wavenumber components of the velocity field due to interpolation. Comparison of this 
spectra to those obtained following the projection steps shows the error resulting from the projec- 
tion schemes. It is also important to note that another factor influencing the errors resulting from 
interpolation and projection is the particle sample size. A smaller ensemble of particles will result 
in a larger error in both the interpolation and projection steps. The simulations presented in §3 
were performed using the same number of particles, 96 3 , as in the tests outlined in Figure 2. 

The constant of proportionality a in (22) and (23) depends on the nature of the particles, i.e., 
actual particles in which each particle in the simulation represents a physical particle, or stochas- 
tic particles in which each particle represents the effect of several. Squires & Eaton (1990) and 
Elghobashi & Truesdell (1993) considered stochastic particles. In that case, it means that for a 
given particle relaxation time and mass loading, the influence of the particles on the fluid motion 
is assumed to be independent of the average distance L between particles. This is fully legitimate 
when L is small with respect to the smallest fluid length scale, the Kolmogorov scale tj. For actual 
particles a = t rd 3 /6, i.e, the volume of a single particle, while for stochastic particles a = a 2 /N p 
where a 2 is the volume fraction of the dispersed phase. For actual particles the volume fraction c* 2 
is determined from the total number of particles in the simulation, a 2 = N p ird 3 / 6/ L\ ox , while for 
stochastic particles, c* 2 is set arbitrarily through the specification of the mass loading <p = a 2 p 2 / Pi- 
It should be noted that the condition L < r] required for correspondence between stochastic and 
actual particles is not met in the previous calculations of two-way coupling by Squires & Eaton 
(1990) and Elghobashi & Truesdell (1993), nor in the present simulations. Thus, the particle re- 
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Table 2: Particle characteristics at zero mass loading. 


(T&stoke, 0069 0-260 0-696 

d/r^ 0-11 0- 21 035 

Re v {4> = 0) 0.38 0.84 1.52 

4(0 = 0) 0.064 0.230 0.580 

4/r fc (0 = 0) 1.26 4.49 11.38 

r e /4(0 = 0) 6.68 1.88 0.74 


laxation time should be considered a parameter of the simulation, along with the number density, 
which is determined based on numerical considerations as shown in Figure 2. The mass loading is 
then changed by varying the material density of the particles. 

2.6 Particle parameters and simulation validation 

The particle parameters are summarized in Table 2. The Stokes relaxation time is denoted (rf 2 ) stokes, 
t f is ^ particle relaxation time obtained from an ensemble and time average over all particles 
with the same material properties. Simulations were performed for three particle relaxation times 
and a series of mass loadings 0 varying from zero (one-way coupling) to unity. Since corrections 
are incorporated for non-linear drag via C D , the particle diameter d and density relative to the 
fluid, pi/ pi (or, equivalently, the relaxation time and density ratio) must be specified. The partic e 
relaxation times were chosen roughly equal to the Kolmogorov, Taylor, and Eulerian integral time 
scales, r* , =0 , r A) , =0 and r^o, of the reference flow at zero mass loading. The corresponding 
diameter variation yields particle Reynolds numbers up to approximately 1.5 for the largest parti- 
cles. For the flow with the intermediate relaxation time, (r p )stokes — 0.2595, at a mass loa mg o 
0 = 0.5, the particles have a diameter d/%= o = 0.156 ^ a relative ^ ensity p2 / pl ~ 3296, tyP ^ Cal 
of gas-solid two-phase flow computations. 

As discussed in §2.3, the fluid velocity locally undisturbed by the particle presence is approx- 
imated in the DNS by the fluid velocity perturbed by all particles. While this is an approximation 
required in order to perform DNS of two-way coupling, it should be pointed out that the inter- 
polation of fluid velocities to particle positions and projection of the coupling force onto the grid 
smooths the influence of the disturbance created by a particle on its own motion. For example, 
using a linear projection to the eight nearest grid points around the particle means that each grid 
node incorporates approximately the velocity disturbance from eight particles. Thus, on average, 


12 



1/8 of the velocity disturbance created by particle n is interpolated back to its position using a 
linear interpolation. For a higher-order accurate interpolation such as that used in this work the 
interpolation of the velocity disturbance created by a particle back to its position is, in general, 
smaller. Increasing the number of particles acts to further reduce the influence of a particle on its 
own motion. 

An estimate of the approximation of the locally undisturbed fluid velocity u", in (10) by U\ t i, 
which is modified by all particles and computed in DNS, was examined through calculation of the 
position and velocity of two groups of particles having the same material properties in the same 
simulation. The following test was undertaken to estimate the effect on second-order statistics from 
approximating the locally undisturbed fluid velocity in (10) by the fluid velocity u Xti which is 
modified by all particles and computed in DNS. The position and velocity of two groups of particles 
having the same material properties were measured in the same simulation. Only the particles of 
the first group influenced the fluid flow. For both groups the drag force was computed with the 
resolved which is, by construction, an approximation of u™ j for the particles of the first group 
but the exact u [ l , for the second group. Statistics from both groups such as the particle kinetic 
energy, the fluid-particle velocity covariance, the fluid kinetic energy along particle trajectories, 
and the fluid-particle energy exchange rate in the fluid kinetic energy equation differed by less 
than one percent. 

3 Results 

The simulations were started from an arbitrary initial condition that was time advanced until the 
rate of energy added to the flow through the forcing balanced the dissipation. Particles were then 
placed randomly throughout the computational domain with an initial velocity identical to the fluid 
velocity at the particle location. Particles were tracked for another development period of roughly 
four relaxation times in order for the particle cloud to reach its own equilibrium condition. Only 
from that point were statistics accumulated by advancing the simulations an additional seven large- 
scale time periods r e . Around six eddy-turnover times were necessary to reach a new equilibrium 
in which production was balanced by both viscous dissipation in the fluid and drag. Statistics were 
then obtained for an additional seven r e . 

3.1 Fluid-phase statistics 

3.1.1 Effect of particles on turbulence statistics 

The equilibrium values of the turbulence kinetic energy, <i ( , and viscous dissipation rate, £\ , in the 
fluid are shown in Figure 3 and 4, respectively. In the absence of gravitational settling small par- 
ticles, d, <C rj, will dissipate turbulence kinetic energy, consistent with the results in Figure 3. The 
results in both Figure 3 and Figure 4 are also in good agreement with Squires & Eaton (1990). Fluid 
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turbulence energy spectra E{k ) are shown in Figure 5. For larger mass loadmg the energy at low 
wavenumbers is diminished independent of t[ 2 , while particles can enhance the high wavenumber 
components, the behaviour at high k depending on <p and r£. This is consistent with the different 
behaviour exhibited by q\ and ^ in Figure 3 and 4. For decreasing rf 2 , small scale motions are 
relatively more energetic compared to the zero mass loading case. For larger r£ and mcreasmg 
<j>, small scale energy reaches a minimum that occurs at lower mass loading that decreases with 
t f Note that for r[ 2 = 0.064 a net production of small scale energy occurs at </>, similar to that 
observed by Squires & Eaton (1990). It is interesting that modulation of the flow by particles is 
capable of modifying the turbulence over the entire spectrum, the degree of modification depend- 
ing on (j) and t[ 2 . The modification is greatest for the smallest particles (t [ 2 = 0.064) which not 
only cause the largest increase in energy at small scales but also affect the lowest wavenumbers. 
These results are contradictory to the notion that particles attenuate, on average, structures having 
a time scale smaller than their relaxation time, which would have lead to a strong damping of high 
wavenumber modes for particles with t [ 2 slightly greater than t*. 

It is also interesting to consider the possibility of “backscatter” of energy from small to large 
scales as a response of the flow to the relative increase in turbulence energy at small scales. In the 
context of this discussion, “backscatter” is simply regarded as an inverse cascade of fluid energy 
at a given wavenumber. Shown in Figure 6 is the evolution with respect to 0 of the fluid transfer 
spectrum T U)1 (fc) for the smallest particles that provide the largest enhancement of small scale 
turbulent motions. The figure shows that, regardless the value of «/», T n ,i{ k ) has a form similar 
to that in single-phase flow with a transfer of energy from large to small scales. The result in 
Figure 6 is similar to that obtained in decaying isotropic turbulence by Elghobashi & Truesdell 
(1993). Finally, it should be pointed out that to completely understand how two-way coupling 
affects non-linear energy transfer in the fluid, a detailed study of triadic interactions is necessary 
(e.g., see Domaradzki et al. 1993). 

The dependence of the Eulerian time macro-scale r £ and Lagrangian integral time scale r* on 
mass loading is shown in Figure 7 and offers further insight into the effect of two-way coupling 
on the different behaviour exhibited by q\ and c t . The Eulerian time macro-scale increaseswith 
larger mass loadmg, eventually reaching a maximum that occurs at larger (j) for increasing t[ 2 . As 
is clear from the figure, the Lagrangian integral time scale exhibits a general increase with (j>. Thus, 
assuming a direct proportionality between these scales, as is often done in turbulence models, is 
not accurate for large mass loadings. 

In decaying turbulence time scales increase with time (e.g., see Chasnov 1994). One can thus 
hypothesise that the more the energy of certain scales decrease, the more their associated time 
scales increase. For loaded flows, Figure 5 shows that small scales are diminished relatively less 
than the large scales. Moreover, r[ characterises more large scale structures than r e , which is more 
representative of the energy cascade. Thus, t € should be augmented less than that t\ which is in 
good accordance with the results. 
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Shown in Figure 8 is the ratio of the particle relaxation time to the Kolmogorov time scale. The 
figure shows that the ratio decreases with increasing mass loading, a consequence of the reduction 
in £\ . This is in turn consistent with the effect of the larger particles on the flow being similar 
to that of the smaller particles as 0 increases. For example, for t( 2 = 0.58 the spectrum E(k ) 
shows a relative increase of small scale energy at high wavenumbers for 0 = 1, similar to that 
observed at smaller loadings for rf 2 = 0.064 (a similar feature may be observed in the spectrum of 
the fluid-particle energy exchange rate shown in Figure 14). The reduction in the time scale ratios 
with increasing 0 complicates interpretation of whether two-way coupling can be more accurately 
described using the ratio of the relaxation time to the time scale of the large eddies, Tj^/r e , or the 
time scale of the smallest eddies, r£/r*. To resolve this issue requires simulations at substantially 
higher Reynolds numbers than can be achieved using DNS in order to provide a much larger 
separation between r e and r*. This may come from an effect induced by the low Reynolds number 
of the simulations that does not improve as the load increases, since diminishes and with it 
T £/ Tjt (see Figure 8). Thus the behaviour of large particles slips towards that of intermediate ones. 

3.1.2 Turbulence transport equations 

For statistically stationary isotropic turbulence modified by momentum exchange with particles, 
the transport equations for the fluid turbulence kinetic energy and dissipation rate are, 

-ei + n 9l + F 91 = 0 , (24) 


£<n - e d2 + n £1 + f Ei — o , 


(25) 


where e dl and e d 2 in (25) are the production by turbulent vortex stretching and viscous destruction 
of dissipation, respectively; e d 2 can be expressed in terms of E(k) as 


poo 

£<12 = 4i^ / k*E{k)dk. (26) 

Jo 

As discussed, for example, in Smith & Reynolds (1991), production by turbulent vortex stretching, 
e dl , is characteristic of a spectral transfer in that it measures the stretching of all turbulent struc- 
tures. The terms F qi and F ei are the contributions from the forcing, n ?1 is the fluid-particle energy 
exchange rate, and n £l the fluid-particle dissipation exchange rate. 


n ?1 — 


<t> 


r 12 


[2 < q\ >2 -912] n £1 = -21/j-p < 

r 12 


0 d(u”^ U 2 ,j) 


dxj 


dxj 


>2 


(27) 


In (27), qi2 is the fluid-particle velocity covariance, < qf >2 is the fluid kinetic energy along the 
particle trajectory, u" t and u 2 i are the fluid and particle velocity fluctuations measured along the 
particle trajectory, respectively. Note that < • > 2 denotes averages over the dispersed phase. 

The production of dissipation £ dl by vortex stretching is shown in Figure 9. The figure shows 
a large reduction with increasing 0 and a weak dependence (reduction) on t[ 2 . The dependence on 
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particle relaxation time appears stronger with increasing 0. Overall, e d i exhibits similar behaviour 
as observed previously with e\ (c.f., Figure 4). The viscous destruction of dissipation e d 2 exhibits 
an interesting behaviour plotted in Figure 10. It initially decreases with increases in loading. For 
the smallest particles, e d2 attains a minimum for ^ > 0.2, which seems to preclude an increase at 
higher (f>. The plateau in e d 2 occurs at higher (j) with increasing t 12- Note also that the plateau in 
e d2 occurs at the same mass loadings as that observed at high wavenumbers in the energy spectra 
in Figure 5. 

Thus the different types of particles do not influence the turbulent quantities in the same way. 
By the relations (26) that link them to the turbulent energy spectrum, one knows that each of them 
measures more specifically certain scales. Thereby this confirms the common expectation that the 
types of particles should play differently on the turbulent scales. 

Shown in Figure 1 1 is the relative behaviour of e dl and e d2> the dissipation of the energy dis- 
sipation rate e £l = £ dl - e d2 , normalized by the quantity it is usually modeled by, -£?/<??• In 
single-phase flows a balance exists between s d i and e d2 and therefore these terms are modeled to- 
gether (e.g., see Smith & Reynolds 1991). The ratio shown in Figure 1 1 increases with (j) except at 
low loadings for large particles where it actually becomes negative, indicating the effect of the par- 
ticles is to provide a source of dissipation. The results in the figure clearly show that modulation of 
the flow by particles can strongly disrupt the equilibrium between e d i and e d2 . This in turn implies 
that £ ei must be modeled differently than in single-phase turbulence, with an explicit dependence 
on (j> and (see Squires & Eaton 1994 for further discussion). 

3.2 Particle source term statistics 

3.2.1 Fluid-particle energy exchange rate 

On average, the particles are an additional dissipation of kinetic energy. The fluid-particle energy 
exchange rate 1 1,, is therefore negative. When normalised by e i,^=o, the evolution of — II 9I dis- 
played in Figure 12 shows that dissipation by drag increases with particle size and mass loading. 
However, the figure also shows that — II 9l /e'i,^=o appears to reach a plateau for larger C>. The in- 
creasingly large contribution of the particles to the total dissipation with larger O is shown more 
clearly in Figure 13 where the ratio of the particle dissipation — II 9l normalised by the total value 
is shown. It should be noted that production of kinetic energy by the forcing, F qi , exhibits a slight 
reduction with increased loading due to the fact that F qi measures the correlation between the 
forcing acceleration and the low wavenumber modes of the velocity which decrease for increasing 
(j>. 

The spectrum of the fluid-particle energy exchange rate, II qi (k), is shown in Figure 14. The 
spectra display two regions, the low-wavenumber portion of the spectrum shows that the fluid tur- 
bulent motion transfers energy to the particles, i.e., the particles act as a sink of kinetic energy. At 
higher wavenumbers the spectrum of the energy exchange rate is positive, indicating that particles 
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are capable of adding kinetic energy to the turbulence. This energy “released” by the particles is 
not immediately dissipated by viscous effects but is in fact responsible for the relative increase of 
small scale energy previously observed in the energy spectra (c.f.. Figure 5a). Thus, in the low 
wavenumber range the fluid drags the particles while at high wavenumbers, particles can drag the 

fluid. 

As shown in Figure 14, the drag of the fluid by the particles is more apparent at smaller relax- 
ation times and for larger loadings, e.g., for r[ 2 = 0.58, shown in Figure 14b, there is essentially no 
wavenumber range over which particles impart kinetic energy to the fluid. The figures also show 
that the absolute value of the negative portion of the spectrum increases weakly with r 12 , indicating 
that the magnitude of If 91 is an increasing function of t[ 2 . This is in accordance with the results 
for U qi shown in Figure 12. For increasing cf>, the region of positive Il 9l ( k ) increases in magnitude 
and the corresponding wavenumber range also shifts towards larger scales. The increase in If qi at 
small scales is not large enough to counter-balance the increase in magnitude at large scales, which 
in turn increases the magnitude of II,, for larger loading, however and as already pointed, falls off 
at higher loads, see Figure 12. 

3.2.2 Fluid-particle dissipation exchange rate 

The source term representing the direct effect of the particles on the dissipation rate, n ei , is shown 
in Figure 15. This quantity undergoes the most striking evolution of the turbulence quantities. For 
small cj) it acts as a sink of dissipation, with the exception of the smallest particles, before becoming 
a source of dissipation as the loading increases. The loading ratio at which If £l changes from sink 
to source also increases with r[ 2 . The evolution of II £l can be more clearly understood from its 
definition: 

/»oo 

n £l = / II ei (k) dk = I 2v 1 k 2 n qi {k) dk. (28) 

Jo J 0 

Particles with small relaxation times drag the fluid at small scales, an effect which increases with 
larger loading. The weighting of the high wavenumbers in (28) accentuates this effect, ultimately 
causing II £I to act as a source of dissipation. 

3.3 Spectral analysis 

The spectral equations for the fluid turbulence, particle fluctuating velocities, and fluid-particle 
covariance are obtained by manipulation of (7) and (10) (e.g., see Baw & Peskin 1971). These 
equations allow one to form the appropriate two-point correlations from which Fourier transfer- 
mation can then be applied to obtain the transport equations in spectral space. 

In homogeneous isotropic turbulence, a projector can be used to pass from directional to tridi- 
mensional quantities. Its application to the spectral equations followed by an integration over 
angular variables yields equations governing the fluid turbulence energy spectrum, Eu(k), the 
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energy spectrum of the fluid turbulence along the particle trajectory, E nfi (k); the fluid-particle co- 
variance spectrum, E l2 {k); and the particle energy spectrum, E 22 {k). For wavenumbers k greater 
than Kp (the radius of the spectral sphere of the forced modes) for which the forcing contribution 
is identically nul these equations are, 


(see 


- 

+ 

+ 

- 

n 9l (fc) = 

Baw & Peskin 1971, Boivin 


T n>1 (k) - 2u 1 k 2 E n {k) + U qi {k) 

(29) 

1 V + ^ ] 

Ei 2 {k) = Ti 2> i(k) + T 12fi (k) 

(30) 

L J 



~~p [ Enpik) -1- <j>E 2 2{k) ] 
T 12 

(31) 

Tooo(k) et [ 2E 22 (k) — Ei 2 (k) ] 

(32) 
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1996). The integrals corresponding to the spectra are. 

(33) 
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(34) 


The terms T n , 1 , T 12 ,i, T 12 , 2 and T 22 , 2 represent non-linear energy transfer in the fluid turbu- 
lence, transfers of the fluid-particle correlated motion by the fluid turbulence along the particle 
path, and by the particle fluctuating motion, and the transfer of particle-particle correlated mo- 
tion by the particle motion, respectively. It should be noted that U qi ( k ) was computed according 
to its formal definition given above in (33). For non-settling particles in homogeneous isotropic 
turbulence, Il 9l (k) can also be expressed in terms of E nfi and E 12 as. 


U qi (k) = -^\E n , 2 (k) - E 12 (k)\ . (35) 

T 12 

The integral of (35) gives (27). In the solution of the equations developed below, E n , 2 {k) is 
approximated by E u (k), corresponding to an assumption of negligible differences in turbulence 
properties measured along particle trajectories compared to those on the grid. For statistically 
stationary flows such as that considered in this work, the time derivatives are zero. Closure of the 
system (29)-(33) then requires models for the transfer terms. 


3.3.1 Closure of the transfer terms - unloaded flows 

For monophase turbulent flow, Pao (1965) proposed a form for T u> i(fc) (typically denoted in 
single-phase turbulence as T(Jfe)) valid in the universal equilibrium range that respects the -5/3 
slope in the inertial sub-range with a stronger damping of energy in the dissipative range. Pao 


18 



(1965) assumed that the energy flux S(k) through wavenumber k is directly proportional to E(k), 
the energy density at wavenumber k. From a dimensional analysis, S(k) can be expressed as, 

S(k)=a~ 1 e\ /3 k 5/3 E(k) with = (36) 

where a is the Kolmogorov constant. For zero mass loading the fluid turbulence spectral equation 
becomes for high Reynolds number turbulence, 

4crV/ 3 fc 5/3 £(/fc) = -2 v x k 2 E{k) , (37) 

dk 

whose solution is, 

E(k) = aei /3 fc- 5/3 exp(-3/2a(A:7 ? ) 4/3 ) . (38) 

Figure 16 compares the theoretical energy and dissipation spectra with those from the simulations 
at zero mass loading. The theoretical result can reproduce quite correctly the simulation results 
provided that a = 3, rather than the usual value of 1.5 valid for high Reynolds number turbulence. 

3.3.2 Closure of the transfer terms - loaded flows 

In two-phase flows with two-way coupling, Baw & Peskin (1971) assumed that I\ i,i. (k) can be 
expressed similarly as in single-phase turbulence in which viscous dissipation is equal to the energy 
transfer rate from the large scales. In an equilibrium forced turbulence the rate of energy transfer 
from the large scales, via the forcing, is dissipated by viscous effects in the fluid or by drag around 
the particles, i.e., F qi = e 1 - n 9l . Thus, the energy transfer rate appropriate for a spectral analysis 
of two-way coupling is the total dissipation, e x - n ?1 , rather than simply £i as in single-phase 
turbulence. Therefore, we propose to replace s j by the energy transfer rate for two-way coupling, 
s l — n 91 , in the expression for S{k). The transfer term Tu .,1 (k) can then be written as, 

T n ,i(k) = - II qi )^k^E u (k ) . (39) 

For Ti 2 ,i(k), a closure analogous to that used in single-phase turbulence is adopted. The terms 
Tiii {k) and T 12 ,i {k) are similar in that T 1U represents the transfer of fluid-fluid correlated motion 
by the fluid turbulence and Ti 2 ,i represents the transfer of fluid-particle correlated motion by the 
fluid. Therefore, it is appropriate to replace E xx in (39) by E v >, while maintaining the same rate of 
energy transfer — II 9l , i.e., 

T n ,i{k) = - n qi )^k 5 ^E l2 (k ) . ( 40 ) 

Baw & Peskin (1971) assumed that, due to particle inertia, Tu,i, 7\ 2) 2 and T 2 2,2 should be very 
small and can be neglected. It seems rather hazardous to imagine particle fluctuating motion similar 
to the fluid turbulence with a cascade of energy, etc. Thus, the development of closure models for 
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T 12>2 and T 22 ) 2, representing non-linear energy transfer by the particle fluctuating motion, would 
require a completely different approach than that typically proposed for T U)1 . Spectral closures for 
Ti 2 2 and T 22i2 are not within the scope of this work and are therefore neglected. 

With the approximations described above, the system of spectral equations becomes 


d\ a-^ei -U qi yE k ^E n (k)] 

k 2 Eu(k) - 

■-^[2E n (k)-En(k)] 

(41) 

dk 
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dla^iei-U^k^Enik)] 


\ + i>\k 2 

En{k) 4 — -pEn(k) 

(42) 

dk 



T 12 


E2i{k) = 

\En(k) 


(43) 


3.3.3 Model without Tn,i 

If one neglects the transfer of the fluid-particle covariance by the fluid turbulence, Ti 2 ,i, then the 
only differential equation to solve is that for E n . The system (41)-(43) in this case is identical to 
that obtained by Baw & Peskin (1971) with - fl 9l instead of e x as the rate of transfer. 


dla-'iex-nj^kWEnik)] 

dk 


= ~2v l k‘ l E ll {k)-^p{2En{k)-E x2 {k)) (44) 

r 12 


E\2{k) — 


1 + k 2 Ti2 


E n (k) 


E22 i,E) — 2 Eu{k) ■ 


(45) 

(46) 


Using (45), the term due to the modulation of the turbulence by particles in the expression (44) for 
Eu can be expressed as 

v ik 2 


4> 


.- [ 2E n (k)-Ei 2 (k)] = -2^ i + ^ 


Yf Eu(k). 


(47) 


The right hand side of (47) is negative definite, indicating that particles act as a sink in (44) across 
the spectrum. This is in contrast, however, to the results presented in §3 which showed that there 
can be an increase in the high wavenumber components of the fluid turbulence, e.g., for the smaller 
particles at higher mass loadings. 


3.3.4 Model with T 12)1 

Inclusion of the closure model for T l2 ,i requires solution of (42). It is also important to point 
out that the term accounting for the modulation of the fluid spectrum, n 9 l (A:), is changed by the 
inclusion of Ti 2>1 . Solving (3 1) for E n and substitution into (35) yields: 


n 


, l(fc ) = - A [ 2 E n (k) - En ik) ] = -20 EuW + 


20 


The contribution in U qi (k) from T n ,i (k) should be positive at high wavenumbers and can therefore 
balance the negative contribution from En(k). This in is in turn consistent with the change of sign 
in U qi (k) observed in Figure 14. Note also that the factor in front of T 12 ,i increases with 0 and 
decreases with r[ 2 . Moreover, the motion of larger particles is less correlated with the fluid, which 
should in turn reduce the fluid-particle covariance transfer (c.f., Eqn. 40). Thus, inclusion of the 
transfer term Ti 2 ,i is appears to behave in accordance with the evolution of the small scale motions 
observed in the DNS. 


3.3.5 Numerical resolution 

The solution of (44)-(46) is more easily accomplished by defining new variables: 

G n (k) = 1/fak 5 ' 3 E n (k), fa = a(e 1 - Yl qi )~ 1/3 (49) 

G u (k) = l/fa fc 5/3 E 12 (k) , fa = «(ei - n 91 )~ 1/3 . (50) 

The general system (41)-(43) becomes for all k greater than Rf, 

dG ^ k l = fa f _ 2 vfa 1 ! 3 - 2-^-/c -5/3 1 G u {k) + fa^k^ 3 G u {k) (51) 

dk I t 12 J r i2 

dGl2(fc i = fo—k-W G u (k) + fa [ ~^iA: 1/3 - -^ fc_5/3 I G 12 W • (52) 

dk r{ 2 L r i2 J 


The system (51)-(52) can be solved numerically for both non-zero and zero values of dG i2 (k)/dk, 
corresponding to inclusion or neglect of T 12) i in the spectral equations. Because En(k) and pre- 
sumably E l2 (k) decrease more rapidly than At 5/3 , conditions on G n {k) and G n (k) are G n (oo) = 
G 12 ( 00 ) - 0. However, these conditions are not convenient for numerical solution. In order to 
facilitate comparison between the DNS results and theoretical predictions, the spectral values at 
k = 4 were used to solve (5 1)-(52). Note from Figure 1 6 that k = 4 corresponds to the beginning 
of the region in which there is good agreement between the DNS results and theory for zero mass 
loading. 

Figure 17 compares the predicted spectra, with and without the closure (40) for T 12>1 , to the 
DNS results for t [ 2 = 0.064 and t [ 2 = 0.58, both at = 1. For t [ 2 = 0.064, inclusion of T n ,\ 
is crucial to obtaining a very good agreement between the DNS results and model prediction. In 
particular, the model prediction does not roll off as rapidly at high wavenumbers and also respects 
the relative increase of small-scale energy observed in the DNS. For t (' 2 = 0.58, the system (44)- 
(46) (neglecting Ti 2> i) is in reasonable agreement with the DNS results and Figure 17 shows that 
inclusion of the model for T 12j1 tends to over estimate E(k). As shown in the figure, however, the 
smaller contribution for rf 2 = 0.58 does not sufficiently damp Ti 2j i for the larger particles. This 
implies that a more accurate model for the transfer of the fluid-particle covariance should have an 
exhibit dependence on r£- 

Figures 18a and 18b demonstrate the effect of T 12]1 on the fluid-particle energy exchange. 
Shown in the figure is the modeled and actual form of the spectrum of the fluid-particle energy 
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exchange rate, U qi {k). The contributions of the terms in (48) are also shown. The figure shows 
that qualitatively that the prediction (48), including the transfer T^i, can reproduce the overall 
behaviour in U qi (k). In particular, the contribution of the transfer of fluid-particle covariance by 
fluid turbulence succeeds in making n gi (7c) positive at high wavenumbers for the small particles. 
The integral value of the positive portion of Tl qi (k) is indeed inversely proportional to t [ 2 and the 
wavenumber at which the effect of the particles changes from a sink to a source of turbulence 
energy also increases with T n » in agreement with DNS results. 

Figure 18 also shows that U qi (k) is over estimated at the lower wavenumbers. One factor 
influencing the over estimation is the neglect of 7 ' 12,2 in the transport equation for E-u- If it is 
assumed that T 12 , 2 is analogous to non-linear transfers in the fluid, T n , 1 and T 12>1 , then T 12>2 would 
be comprised of two portions, negative at low wavenumbers and positive at higher wavenumbers. 
In the spectral equation (31) for E n {k), T 12)2 would then be a sink at low wavenumbers and a 
source at high wavenumbers. At the large scales this would result in a reduction in Eyiik) and, 
according to (35), a decrease of n gi (A:) at low wavenumbers. Similarly, the relative mcrease in 
II 9l ( k ) at high wavenumbers due to the contribution from Ti 2j2 would also reduce the discrepancy 
between model predictions and DNS results in Figure 18. 

4 Conclusion 

Direct numerical simulation of the incompressible Navier-Stokes equations has been used to in- 
vestigate turbulence modulation by particles in isotropic turbulence within micro-gravity environ- 
ments. The flow field was forced at the low wavenumbers to maintain a statistically stationary 
condition. Three relaxation times ranging from the Kolmogorov time scale to the Eulenan time 
scale of the reference fluid flow (<p = 0) and loading ratios </> ranging from 0 to 1 comprised the 
particle parameter space. 

For non-zero mass loading, particles increasingly dissipate fluid kinetic energy as the loading 
ratio increases, with the reduction in the kinetic energy being relatively independent of the relax- 
ation time. Simultaneously, viscous dissipation in the fluid decreases with increases in <f), being 
larger for particles with smaller relaxation times. Furthermore, the ratio of the kinetic energy to the 
dissipation rate, that defines the Eulerian time macro-scale t e , differs noticeably from that of the 
fluid Lagrangian integral time scale t(, which increases with <p. The different response of quan- 
tities such as the kinetic energy and dissipation rate to increased loading and changes in particle 
relaxation time is in turn linked to the fluid turbulence energy spectra. DNS results show that 
there is a non-uniform distortion of E(k), with a relative increase in small-scale energy. The non- 
uniform distortion significantly affects the transport of E\ since the production of dissipation e^i 
and destruction of dissipation E& exhibit completely different behaviours. For example, for small 
relaxation times and large mass loadings, particles can be a source of dissipation, rather than a sink 
as conventionally modeled. 
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The fluid-particle energy exchange rate, n^, increases relative to the total dissipation for both 
larger loading and relaxation time. The spectrum H qi (k) shows that, while the fluid drags the 
particles at low wavenumbers, particles with small r ^ drag the fluid increasingly at larger loading. 
Thus, particles “release” to turbulent small scale motions part of the energy extracted from the 
larger scales. This phenomenon in turn explains the switch from sink to source of the fluid-particle 
dissipation exchange rate U £l that occurs at lower (j) for small particles. A spectral analysis shows 
that the increase of the high wavenumber portion of the fluid energy spectrum can be attributed to 
transfer of the fluid-particle covariance by the fluid turbulence. 

While the approach used in the present study— DNS— is appropriate for detailed analyses of 
particle-turbulence interactions, there are issues relevant to two-way coupling which cannot be re- 
solved using direct simulations. The low Reynolds numbers and limited range of scales in DNS, 
for example, prevent a determination of whether two-way coupling is best described in terms of 
large- or small-scale variables, i.e., in terms of T^/r e or Because of the “global” distor- 

tion of the turbulence across the entire spectrum, a description of two-way coupling in terms of 
small-scale variables may not be the most appropriate, as is conventionally assumed. Calculations 
at substantially higher Reynolds numbers are required in order to obtain a wide separation between 
the energy-containing and dissipating scales. While DNS of single-phase (homogeneous) turbu- 
lence can be performed at higher Reynolds numbers, the computational constraints, e.g., adequate 
particle sample sizes, necessary for accurate resolution of two-way coupling virtually prohibits the 
use of DNS. 


23 


References 

[1] Balachandar, S. & Maxey, M.R. 1989 Methods for evaluating fluid velocities in spectral sim- 
ulations of turbulence. J. Comp. Physics, vol. 83, pp. 96. 

[2] Baw, P.S.H. & Peskin, R.L. 1971 Some aspects of Gas-Solid Suspension Turbulence. J. of 
Basic Engineering, pp. 63 1 . 

[3] Boivin, M. 1996 Etude de l’lnfluence des Particules sur la Turbulence a partir de Simulations 
Directes et de Simulations des Grandes Echelles d’ecoulements Diphasiques Gaz-Solides 
Homogenes Isotropes Stationnaires, Ph.D. thesis. 

[4] Clift, R., Grace, J.R. & Weber, M.E. 1978 Bubbles, Drops, and Particles. Academic Press, 

pp. 106. 

[5] Comte-Bellot, G. & Corrsin, S. 1971 Simple EulerianTime Correlation of Full-and Narrow- 
Band Velocity Signals in Grid-Generated Isotropic Turbulence. J. Fluid Meek, vol. 186, pp. 
199. 

[6] Crowe, C.T. 1982 Review-numerical methods for dilute gas-particle flows. J. Fluids Engng., 
vol. 104, pp. 297. 

[7] Crowe, C.T., Trott, T.R. & Chung, J.N. 1 996 Numerical models for two-phase turbulent flows. 
Ann. Rev. Fluid Mech ., vol. 28, pp. 1 1 . 

[8] Csanady, G.T. 1963 Turbulent diffusion of heavy particles in the atmosphere. J. Atm. Sci ., 
vol. 20, pp. 201. 

[9] Deutsch, E. & Simonin, O. 1991 Large Eddy Simulation applied to the modeling of partic- 
ulate transport coefficients in turbulent two-phase flows. Proc. 8th Int. Symp. on Turbulent 
Shear Flows, (Univ. of Munich), vol. 1, pp. 1011. 

[10] Domaradzki, J., Liu, W. & Brachet, M. 1993 An analysis of subgrid-scale interactions in 
numerically simulated isotropic turbulence. Phys. Fluids A, vol. 5, pp. 1747. 

[1 1] Elghobashi, S.E. & Truesdell, G.C. 1993 On the two-way Interaction between Homogeneous 
Turbulence and Dispersed Solid Particles. I: Turbulence Modification. Phys. Fluids A, vol. 5, 
pp. 1790. 

[12] Eswaran, V. & Pope, S.B. 1988 An Examination of Forcing in Direct Numerical Simulations 
of Turbulence. Computers & Fluids, vol. 16, No 3, pp. 257. 

[13] Fleckhaus D„ Hishida, K. & Maeda, M. 1987 Effect of laden solid particles on the turbulent 
flow structure of a round free jet. Exp. Fluids, vol. 5, no. 5, pp. 323. 


24 


[14] Gatignol, R. 1983 The Faxen formulae for a rigid particle in an unsteady non-uniform Stokes 
flow. J. de Mec. Th. et Appl., vol. 1, pp. 143. 

[15] Gore, R.A. & Crowe, C.T. 1989 Effect of Particle Size on Turbulence Intensity. Int. J. Multi- 
phase Flow, vol. 15 (2), pp. 279. 

[16] Hardalupas, Y., Taylor A.M.K.P. & Whitelaw, J.H. 1989 Velocity and particle-flux charac- 
teristics of turbulent particle-laden jets. Proc. Roy. Soc., Series A, vol. 426, no. 1870, pp. 
31. 

[17] Hinze, J.O. 1975 Turbulence McGraw-Hill Book Company, New-York. 

[18] Koch, D. 1990 Kinetic theory for a monodisperse gas-solid suspension. Phys. Fluids, vol. 
2(1), pp. 1711-1723. 

[19] Kulick, J.D., Fessler, J.R. & Eaton J.K. 1994 Particle response and turbulence modification 
in fully developed channel flow. J. Fluid Mech., vol. 277 , pp. 109. 

[20] Maxey, R. & Riley, J. 1983 Equation of Motion for a Small Rigid Sphere in a Turbulent Fluid 
Flow. Phys. Fluids, vol. 26, pp. 883. 

[21] Modarress, D., Tan, H. & Elghobashi, S.E. 1984 Two component LDA measurement in a 
two-phase turbulent jet. AIAA J., vol. 22, no. 5, pp. 624. 

[22] Pao, Y.H. 1965 Structure of turbulent velocity and scalar fields at large wavenumbers, phase 
fluide , vol. 8, pp. 1063. 

[23] Rogallo, R.S. 1981 NASA Technical Memo 81315. 

[24] Rogers, C.B. & Eaton, J.K. 1991 The effect of small particles on fluid turbulence in a flat- 
plate turbulent boundary layer in air. Phys. Fluids, vol. 3, no. 7, pp. 928. 

[25] Saffman, P.G. 1973 On the settling speed of free and fixed suspension. Stud. Appl. Math., vol. 
LII(2), pp. 115-127. 

[26] Shuen, J-S., Solomon, A.S.P., Zhang, Q-F. & Faeth, G.M. 1985 Structure of Particle-Laden 
Jets: Measurements and Predictions. AIAA J., vol. 23, no. pp. 396. 

[27] Simonin, O., Deutsch, E. & Minier, J.P. 1993 Eulerian Prediction of the Fluid/Particle Corre- 
lated Motion in Turbulent Two-Phase Flows. Applied Scientific Research, vol. 5, pp. 275. 

[28] Simonin, O., Deutsch, E., Boivin, M., 1995, Large Eddy Simulation and Second-Moment 
Closure Model of Particle Fluctuating Motion in Two-Phase Turbulent Shear Flows. In Se- 
lected Papers from the Ninth Int. Symp. on Turbulent Shear Flows, F. Durst, N. Kasagi, B.E. 
Launder, F.W. Schmidt, K. Suzuki, J.H. Whitelaw (Editors), Springer- Verlag, pp. 85. 


25 



[29] Smith, L.M. & Reynolds, W.C. 1991 The Dissipation-Range Spectrum and the Velocity 
Derivative Skewness in Turbulent Flows. Phys. Fluids A, vol. 3, pp. 992. 

[30] Squires, K.D. & Eaton, J.K. 1990 Particle response and turbulence modification in isotropic 
turbulence. Phys. Fluids A, vol. 2, p. 1191. 

[31] Squires, K.D. & Eaton, J.K. 1991a Preferential concentration of particles by turbulence. 
Phys. Fluids A, vol. 3, pp. 1 169. 

[32] Squires, K.D. & Eaton, J.K. 1991b Measurements of particle dispersion obtained from direct 
numerical simulations of isotropic turbulence. J. Fluid Mech., vol. 226, pp. 1 . 

[33] Squires, K.D. & Eaton, J.K. 1994 Effect of selective modification of turbulence on two- 
equation models for particle-laden turbulent flows. Journal of Fluids Engineering, vol. 1 16 
(4),pp. 778. 

[34] Unverdi, S.O. & Tryggvason, G. 1992 A front-tracking method for viscous, incompressible, 
multi-fluid flows. J. Comp. Physics, vol. 100(1), pp. 25. 

[35] Tsuji, Y., Morikawa, Y. & Shiomi, H. 1984 LDV measurements of an air-solid two-phase 
flow in a vertical pipe. J. Fluid Mech., vol. 139, pp. 417. 

[36] Wang, L.P. & Maxey, M.R. 1993 Settling Velocity and Concentration Distribution of Heavy 
Particles in Homogeneous Isotropic Turbulence. J. Fluid Mech., vol. 256, pp. 27. 

[37] Yeung, P.K. & Pope, S.B. 1988 An algorithm for tracking fluid particles in numerical simu- 
lations of homogeneous turbulence. J. Comp. Physics , vol. 79, pp. 373. 

[38] Yeung, P.K. & Pope, S.B. 1989 Lagrangian Statistics from Direct Numerical Simulations of 
Isotropic Turbulence. J. Fluid Mech., vol. 207, pp. 531. 


26 




Figure 2: Turbulence energy spectrum from one realisation of a statistically stationary flow at zero 

mass loading. spectrum before interpolation; after interpolation on to a mesh shifted 

half a grid cell; after interpolation on to 96 3 random particle locations and projection back on 

to the grid using the PIC method (22); after interpolation on to 96 3 random particle locations 

and projection back on to the grid using linear projection (23). 
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Figure 12: Particle dissipation normalised by fluid viscous dissipation U qi M,*=o 
0.064 ; Tj2 = 0.23 ; A-A = 0.58. 
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Figure 13: Particle dissipation normalised by total dissipation -U q J(ei~ II 9l ) . c 

0 — □ r£ = 0.23 ; A -A rg = 0.58. 


















